Análisis exploratorio de datos

Vamos a hacer EDA (Exploratory Data Analysis) de datos municipios de Brasil de 2019 de tres fuentes:

• Pagos del plan Bolsa Familia del Portal de Transparencia
• El Instituto Brasileño de Geografía y Estadística (IBGE)
• El Instituto de Investigación Económica Aplicada

Vamos a suponer que nuestro objetivo final es predecir el PIB a nivel municipio. Esto podría ser útil, por ejemplo, si sabemos que para algunos municipios este dato es poco confiable, mientras para otros no. También si es un dato que se publica con rezago con respecto a otros que están disponibles más rápido.

WARNING: esto es solamente un escenario hipotético. Si de verdad quisiéramos hacer este ejercicio deberíamos analizar muchos aspectos, como:

• ¿Cómo se estima el PIB? ¿Hay algún método más simple que me permita estimar el dato sin machine learning?
• ¿La población de municipios sin dato confiable tiene una distribución similar en sus features a la población con dato confiable? En tal caso, ¿podemos suponer que el DGP (proceso generador de datos) es el mismo?
• ¿Los features van a estar disponibles al momento de la predicción?
• ¿Hay otros features importantes disponibles?
• Etc.

WARNING 2: en este análisis nuestro objetivo es explorar y entender los datos, formular hipótesis (por ejemplo, acerca de qué features pueden ser importantes), etc.. Esto quiere decir las transformaciones que hagamos acá no necesariamente vayan a ser las mismas que apliquemos cuando entrenemos un modelo.

WARNING 3: el siguiente es un EDA de entre muchos posibles. ¡Seguramente no sea el mejor porque es para dar clase!

library(tidyverse)
library(janitor)
library(skimr)
library(GGally)
library(ggpubr)
library(hexbin)
library(glue)
library(scales)
# library(kableExtra) # dont load it because it breaks skimr

# install.packages("remotes")
# remotes::install_github("rpkgs/gg.layers")
library(gg.layers)

Levantamos los tres datasets:

file_bolsa = "data/working/bolsafamilia_municipios_2019.csv"
file_ibge = "data/working/ibge_municipios_2019.csv"
file_varios = "data/working/varios_municipios_2019.csv"

df_bolsa = read_csv(file_bolsa, show_col_types=F)
df_ibge = read_csv(file_ibge, show_col_types=F)
df_varios = read_csv(file_varios, show_col_types=F)

Inspección inicial

Inspeccionamos algunas características básicas de los data.frames:

glimpse(df_bolsa)

## Rows: 5,570
## Columns: 4
## $ uf                     <chr> "GO", "MG", "GO", "MG", "PA", "CE", "BA", "BA",~
## $ num_beneficiarios_mean <dbl> 995.2500, 273.6667, 685.8333, 539.0000, 34341.5~
## $ suma_valor             <dbl> 1973921, 461869, 1164947, 866639, 84271320, 513~
## $ municipio              <chr> "abadia_de_goias", "abadia_dos_dourados", "abad~

head(df_bolsa)

## # A tibble: 6 x 4
##   uf    num_beneficiarios_mean suma_valor municipio          
##   <chr>                  <dbl>      <dbl> <chr>              
## 1 GO                      995.    1973921 abadia_de_goias    
## 2 MG                      274.     461869 abadia_dos_dourados
## 3 GO                      686.    1164947 abadiania          
## 4 MG                      539      866639 abaete             
## 5 PA                    34342.   84271320 abaetetuba         
## 6 CE                     1852.    5131327 abaiara

• UF (Unidad Federativa) indica el estado.
• num_beneficiarios_mean es el número de beneficiarios promedio a lo largo del año.
• suma_valor es la suma de los pagos totales en el año.

glimpse(df_ibge)

## Rows: 5,570
## Columns: 14
## $ nome_da_grande_regiao                   <chr> "Centro-oeste", "Sudeste", "Ce~
## $ codigo_uf                               <dbl> 52, 31, 52, 31, 15, 23, 29, 29~
## $ nome_da_unidade_da_federacao            <chr> "Goiás", "Minas Gerais", "Goiá~
## $ codigo_municipio                        <dbl> 5200050, 3100104, 5200100, 310~
## $ semiarido                               <chr> "Não", "Não", "Não", "Não", "N~
## $ hierarquia_urbana_principais_categorias <chr> "Metrópole", "Centro Local", "~
## $ vab_agro                                <dbl> 5536.665, 32756.932, 56244.012~
## $ vab_industria                           <dbl> 37221.827, 12233.404, 24158.34~
## $ pib                                     <dbl> 202315.97, 127581.83, 346767.5~
## $ pib_percapita                           <dbl> 23061.21, 18254.66, 17302.04, ~
## $ atividade_maior_vab                     <chr> "Demais serviços", "Administra~
## $ grado_urbanizacion                      <dbl> 96.93234, 51.44562, 63.75355, ~
## $ tipo_urbano                             <chr> "Urbano", "RuralAdjacente", "R~
## $ municipio                               <chr> "abadia_de_goias", "abadia_dos~

glimpse(df_varios)

## Rows: 5,596
## Columns: 9
## $ area_geografica_2019                  <dbl> 147.256, 881.064, 1045.127, 1817~
## $ despesa_ciencia_tecnologia_2019       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ despesa_educacao_cultura_2019         <dbl> 13291163, 6061623, 11262623, 124~
## $ despesa_pessoal_encargos_sociais_2019 <dbl> 18676874, 11345153, 25495707, 26~
## $ despesa_transporte_2019               <dbl> 2192450.55, 1294621.67, 946773.6~
## $ exportacoes_2019                      <dbl> 538966, NA, NA, NA, 70309528, NA~
## $ importacoes_2019                      <dbl> 8894, NA, NA, NA, 42478, NA, NA,~
## $ populacao_2019                        <dbl> 8773, 6989, 20042, 23237, 157698~
## $ municipio                             <chr> "abadia_de_goias", "abadia_dos_d~

• despesa_pessoal_encargos_sociais indica gasto en personal y cargas sociales.
• el nombre del municipio municipio ya está normalizado entre los datasets.

Chequeo de duplicados y joins

Vamos a hacer un join de los datasets usando municipio como clave (ya sabemos que el campo está normalizado). Antes de hacerlo es conveniente revisar si hay observaciones duplicadas o si algún dataset tiene observaciones (municipios) faltantes.

dups_bolsa = df_bolsa %>% get_dupes("municipio")

## No duplicate combinations found of: municipio

dups_ibge = df_ibge %>% get_dupes("municipio")

## No duplicate combinations found of: municipio

dups_varios = df_varios %>% get_dupes("municipio")

## No duplicate combinations found of: municipio

munis_bolsa = df_bolsa %>% pull(municipio) %>% unique()
munis_ibge = df_ibge %>% pull(municipio) %>% unique()
munis_varios = df_varios %>% pull(municipio) %>% unique()

cat(
  length(munis_bolsa),
  length(munis_ibge),
  length(munis_varios)
)

## 5570 5570 5596

¿Los de df_bolsa y df_ibge son los mismos?

# intersect(munis_bolsa, munis_ibge)
length(intersect(munis_bolsa, munis_ibge))

## [1] 5570

¿Cuáles de df_varios no están en df_bolsa/df_ibge?

setdiff(munis_varios, munis_ibge)

##  [1] "assunguy_de_cima"              "augusto_severo"               
##  [3] "barcellos"                     "cimbres"                      
##  [5] "cococi"                        "conchas_2"                    
##  [7] "entre_rios_3"                  "mecejana"                     
##  [9] "monsaras"                      "muribeca_2"                   
## [11] "nossa_senhora_do_o_de_goyanna" "olivenca_2"                   
## [13] "palmira"                       "ponte_do_itabapoana"          
## [15] "porto_de_cima"                 "riacho"                       
## [17] "santa_cristina_do_pinhal"      "santa_izabel"                 
## [19] "santo_amaro_2"                 "sao_goncalo_2"                
## [21] "sao_joao_marcos"               "sao_miguel_2"                 
## [23] "sao_sebastiao_do_parahyba"     "trancoso"                     
## [25] "vertentes_2"                   "vila_franca"                  
## [27] "villa_verde"

¿Cuáles de df_bolsa/df_ibge no están en df_varios?

setdiff(munis_ibge, munis_varios)

## [1] "campo_grande_3"

¿Cuáles están en todos los dfs?

# aplica la funcion intersect recursivamente sobre la lista
munis_en_comun = Reduce(
  intersect, list(munis_bolsa, munis_ibge, munis_varios))
length(munis_en_comun)

## [1] 5569

# Reduce(fn, elementos)

---

Hacemos el join. Es importante saber qué observaciones estamos conservando o eliminando en cada paso.

df = df_bolsa %>%
  full_join(df_ibge, by="municipio") %>% # 5570 municipios 
  left_join(df_varios, by="municipio") 
  # excluimos los que estan en varios que no estan en los otros dos (27 municipios)

Renombrar variables

Renombramos algunas columnas por comodidad:

df = df %>%
  rename(
    nome_regiao = nome_da_grande_regiao
    ,nome_uf = nome_da_unidade_da_federacao
    ,hierarquia = hierarquia_urbana_principais_categorias
    ,tipologia = tipo_urbano
  )

Para normalizar o limpiar los nombres de las columnas es muy útil la función janitor::clean_names() (en este caso ya la aplicamos cuando generamos los datasets).

Esta librería también tiene otras funciones muy útiles, como remove_empty() para eliminar filas o columnas vacías o remove_constant() para eliminar filas o columnas constantes.

Vamos a eliminar el sufijo del año en los nombres de las columnas porque sabemos que todos los datos son del mismo año:

df = df %>% rename_all(function(x) str_remove(x, "_2019"))

Un vistazo

La librería skimr nos da un vistazo general de los datos.

skim(df)

Data summary

Name	df
Number of rows	5570
Number of columns	25
_______________________
Column type frequency:
character	8
numeric	17
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
uf	0	1	2	2	0	27	0
municipio	0	1	3	32	0	5570	0
nome_regiao	0	1	3	12	0	5	0
nome_uf	0	1	4	19	0	27	0
semiarido	0	1	3	3	0	2	0
hierarquia	0	1	9	18	0	5	0
atividade_maior_vab	0	1	10	84	0	10	0
tipologia	5	1	6	22	0	5	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
num_beneficiarios_mean	0	1.00	2474.52	9208.35	2.00	346.71	947.50	2439.04	4.513838e+05	▇▁▁▁▁
suma_valor	0	1.00	5594117.72	18428002.57	3216.00	704285.50	2101369.00	5734926.75	8.539214e+08	▇▁▁▁▁
codigo_uf	0	1.00	32.38	9.83	11.00	25.00	31.00	41.00	5.300000e+01	▂▅▇▅▂
codigo_municipio	0	1.00	3253590.77	984910.34	1100015.00	2512125.75	3146280.00	4119189.50	5.300108e+06	▂▅▇▅▂
vab_agro	0	1.00	55783.48	101866.43	0.00	10113.36	25638.79	59948.02	1.575333e+06	▇▁▁▁▁
vab_industria	0	1.00	248797.85	1429430.02	449.30	4294.51	13302.82	72493.60	5.735987e+07	▇▁▁▁▁
pib	0	1.00	1326594.43	12696536.12	15515.27	87586.45	190970.75	510100.68	7.638060e+08	▇▁▁▁▁
pib_percapita	0	1.00	24546.81	25552.22	4482.85	10515.57	18182.28	29882.71	4.648835e+05	▇▁▁▁▁
grado_urbanizacion	5	1.00	48.28	34.28	0.00	0.00	54.78	77.68	9.967000e+01	▇▂▆▆▆
area_geografica	12	1.00	1526.55	5610.80	3.56	204.90	417.84	1025.39	1.595333e+05	▇▁▁▁▁
despesa_ciencia_tecnologia	289	0.95	63972.80	1805472.41	0.00	0.00	0.00	0.00	1.188419e+08	▇▁▁▁▁
despesa_educacao_cultura	289	0.95	33200752.21	214626480.31	0.00	5704397.12	11392275.52	24730333.47	1.358956e+10	▇▁▁▁▁
despesa_pessoal_encargos_sociais	286	0.95	64558552.22	463388560.08	278944.70	10610257.70	18681271.24	39789851.65	2.508161e+10	▇▁▁▁▁
despesa_transporte	289	0.95	2866400.07	72043623.38	0.00	5720.00	483191.64	1500556.47	5.171955e+09	▇▁▁▁▁
exportacoes	3379	0.39	102088499.29	484266387.22	1.00	350150.50	4824744.00	41587247.00	1.116350e+10	▇▁▁▁▁
importacoes	3553	0.36	87857804.51	514210623.13	48.00	84013.00	1237977.00	12267004.00	1.012875e+10	▇▁▁▁▁
populacao	13	1.00	37600.16	221267.52	781.00	5446.00	11638.00	25501.00	1.225202e+07	▇▁▁▁▁

skim_eda = partition(skim(df))
names(skim_eda)

## [1] "character" "numeric"

Pueden revisar la vignette de skimr para conocer más opciones.

Análisis de datos faltantes

Calculamos el % de datos faltantes por variable:

map_dbl(df, function(x) mean(is.na(x)) * 100) %>% 
  sort(decreasing=T) %>% 
  as.data.frame()

##                                            .
## importacoes                      63.78815081
## exportacoes                      60.66427289
## despesa_ciencia_tecnologia        5.18850987
## despesa_educacao_cultura          5.18850987
## despesa_transporte                5.18850987
## despesa_pessoal_encargos_sociais  5.13464991
## populacao                         0.23339318
## area_geografica                   0.21543986
## grado_urbanizacion                0.08976661
## tipologia                         0.08976661
## uf                                0.00000000
## num_beneficiarios_mean            0.00000000
## suma_valor                        0.00000000
## municipio                         0.00000000
## nome_regiao                       0.00000000
## codigo_uf                         0.00000000
## nome_uf                           0.00000000
## codigo_municipio                  0.00000000
## semiarido                         0.00000000
## hierarquia                        0.00000000
## vab_agro                          0.00000000
## vab_industria                     0.00000000
## pib                               0.00000000
## pib_percapita                     0.00000000
## atividade_maior_vab               0.00000000

y por municipio:

na_by_row = df %>% apply(1, function(x) mean(is.na(x)) * 100)
names(na_by_row) = df$municipio
na_by_row %>% sort(decreasing=T) %>% head(10)

##  balneario_rincao  mojui_dos_campos paraiso_das_aguas    pescaria_brava 
##                40                40                40                40 
##    pinto_bandeira    campo_grande_3           conchas      entre_rios_2 
##                40                32                32                32 
##          muribeca           nazaria 
##                32                32

mean(na_by_row > 0)

## [1] 0.7204668

También podemos recurrir a skimr:

skim_eda$numeric %>% 
  select(skim_variable, n_missing) %>% 
  filter(n_missing > 0)

Variable type: numeric

skim_variable	n_missing
grado_urbanizacion	5
area_geografica	12
despesa_ciencia_tecnologia	289
despesa_educacao_cultura	289
despesa_pessoal_encargos_sociais	286
despesa_transporte	289
exportacoes	3379
importacoes	3553
populacao	13

skim_eda$character %>% 
  select(skim_variable, n_missing) %>% 
  filter(n_missing > 0)

Variable type: character

skim_variable	n_missing
tipologia	5

En un caso real hay que ser muy cuidadoso con los datos faltantes. El tratamiento que les demos depende del tipo de modelo predictivo que usemos: en algunos casos vamos a necesitar imputarlos mientras que ciertas metodologías admiten NAs.

Durante el EDA tomaremos las decisiones que más nos sirvan para entender los datos y comunicar los resultados. Para esto es fundamental preguntarse por la interpretación de los missing variable por variable.

En este ejemplo vamos a suponer que los NA en las variables monetarias indican que el monto es igual a 0.

vars_monetarias = c(
  "despesa_ciencia_tecnologia", "despesa_educacao_cultura", 
  "despesa_pessoal_encargos_sociais", "despesa_transporte",
  "exportacoes", "importacoes", "pib", "pib_percapita", 
  "suma_valor", "vab_agro", "vab_industria"
)
df = df %>% 
  mutate_at(vars_monetarias, function(x) ifelse(is.na(x), 0, x))

Además vemos que los datos de población con NA se pueden reconstruir con el PIB y el PIB per cápita.

df = df %>% 
  mutate(
    populacao = ifelse(
      is.na(populacao), pib / pib_percapita * 1000, populacao)
    ,populacao = round(populacao, 0)
  )

Creación de variables

Para explorar los datos puede ser útil crear nuevas columnas o transformar las que tenemos.

df = df %>%
  mutate(
    densidad = populacao / area_geografica
    ,suma_perbenef = suma_valor / num_beneficiarios_mean
  ) 

df = df %>% 
    mutate_at(
    vars(
      vab_agro, vab_industria, despesa_ciencia_tecnologia
      ,despesa_pessoal_encargos_sociais, despesa_educacao_cultura
      ,despesa_transporte, exportacoes, importacoes
      ,num_beneficiarios_mean, suma_valor
    ), list(pc=function(x) x / .$populacao)
  )

# recodificar variables
df = df %>% 
  mutate_at(vars(codigo_municipio, codigo_uf), as.character)

# categorías nuevas en función de valores de otras variables
df = df %>% 
  mutate(
    semiarido = ifelse(semiarido == "Sim", T, F),
    regiao = case_when(
      nome_regiao == "Norte" ~ "N"
      ,nome_regiao == "Nordeste" ~ "NE"
      ,nome_regiao == "Centro-oeste" ~ "CO"
      ,nome_regiao == "Sudeste" ~ "SE"
      ,nome_regiao == "Sul" ~ "S"
      ,TRUE ~ "error"
    )
  )

Podemos usar forcats::fct_lump() para agrupar categorías poco frecuentes.

df$atividade_maior_vab %>% unique()

##  [1] "Demais serviços"                                                                     
##  [2] "Administração, defesa, educação e saúde públicas e seguridade social"                
##  [3] "Agricultura, inclusive apoio à agricultura e a pós colheita"                         
##  [4] "Indústrias de transformação"                                                         
##  [5] "Pecuária, inclusive apoio à pecuária"                                                
##  [6] "Eletricidade e gás, água, esgoto, atividades de gestão de resíduos e descontaminação"
##  [7] "Comércio e reparação de veículos automotores e motocicletas"                         
##  [8] "Produção florestal, pesca e aquicultura"                                             
##  [9] "Indústrias extrativas"                                                               
## [10] "Construção"

df %>% 
  mutate(
    atividade = forcats::fct_lump(
      atividade_maior_vab, prop=0.1, other_level="Otros")
  ) %>% 
  pull(atividade) %>% 
  unique()

## [1] Demais serviços                                                     
## [2] Administração, defesa, educação e saúde públicas e seguridade social
## [3] Agricultura, inclusive apoio à agricultura e a pós colheita         
## [4] Otros                                                               
## 4 Levels: Administração, defesa, educação e saúde públicas e seguridade social ...

---

A veces es útil separar el dataset según el tipo de variables una vez que terminamos el preprocesamiento.

df_num = df %>% select_if(is.numeric)
df_cat = df %>% select_if(function(x) !is.numeric(x))

Análisis de variables numéricas

Usamos histogramas y densidades para visualizar la distribución.

ggplot(df, aes(x=pib_percapita)) +
  geom_histogram(alpha=0.5) +
  NULL

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(df, aes(x=log10(pib_percapita))) +
  geom_histogram(alpha=0.5) +
  NULL

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

La escala logarítmica es una escala no lineal en la que moverse una unidad implica multiplicarse por las base del logaritmo (2, 10 o e, por ejemplo). \(\text{log}_{10}\) mide, por ejemplo, órdenes de magnitud (\(x\) difiere en un orden de magnitud de \(y\) si \(x\) es aproximadamente 10 veces mayor que \(y\)).

Entonces, la transformación logarítmica es útil cuando la variabilidad relevante está en el orden de magnitud, no en la escala original. Esto es común en el caso de variables monetarias, cuando las variaciones relevantes son porcentuales y no nominales.

Otra manera de visualizar lo mismo:

ggplot(df, aes(x=pib_percapita)) +
  geom_histogram(alpha=0.5) +
  scale_x_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format("log10", math_format(10^.x))
  ) +
  NULL

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Podemos comparar la distribución de una variable entre grupos.

ggplot(df, aes(x=log(pib_percapita), fill=semiarido)) +
  geom_histogram(alpha=0.5) +
  NULL

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# cuando usamos density normalizamos la distribucion en cada grupo
ggplot(df, aes(x=log(pib_percapita), fill=semiarido)) +
  geom_density(alpha=0.5, adjust=2) +
  NULL

# adjust determina el nivel de suavidad
# adjust=2 --> bandas de 2 veces el ancho default

Con la librería ggridges podemos hacer algunos gráficos lindos:

ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill =..x..)) +
  ggridges::geom_density_ridges_gradient(scale=1) +
  scale_fill_viridis_c() +
  NULL

## Picking joint bandwidth of 0.188

ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill=tipologia)) +
  ggridges::stat_density_ridges(
    quantile_lines=T, quantiles=c(.25,.5,.75), alpha=0.7) +
  NULL

## Picking joint bandwidth of 0.188

# a veces es útil visualizar la cantidad de obs al mismo tiempo!
ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill=tipologia)) +
  ggridges::geom_density_ridges(jittered_points=T, point_size=0.5) +
  NULL

## Picking joint bandwidth of 0.188

También podemos usar boxplots:

Fuente: https://r4ds.had.co.nz/

ggplot(df, aes(
  x=fct_reorder(nome_regiao, pib_percapita, .fun="median"), 
  y=pib_percapita, fill=nome_regiao)) +
  geom_boxplot() + 
  scale_y_continuous(trans='log10') +
  # facet_wrap(~ tipologia) +
  theme_minimal() +
  NULL

# eliminando los outliers, invirtiendo las coordenadas, y mostrando los puntos:
ggplot(df, aes(
  x=fct_reorder(hierarquia, pib_percapita, .fun="median"), 
  y=pib_percapita, fill=hierarquia)
) +
  geom_jitter(aes(color=hierarquia), alpha=0.2, size=0.5) +
  gg.layers::geom_boxplot2(alpha=0.5) + 
  scale_y_continuous(trans='log10') +
  # facet_wrap(vars(nome_regiao)) +
  coord_flip() +
  theme_minimal() +
  theme(legend.position="none") +
  NULL

Podemos calcular estadísticos globales o por grupo:

quantile(df_num$pib_percapita, seq(0,1,0.1))

##         0%        10%        20%        30%        40%        50%        60% 
##   4482.850   8124.319   9620.964  11579.717  14544.326  18182.285  22222.446 
##        70%        80%        90%       100% 
##  26959.432  33584.336  45617.949 464883.490

# para todas las variables:
list_cuantiles = df_num %>% 
  map(function(x) quantile(x, seq(0, 1, 0.25), na.rm=T) )

list_cuantiles %>% 
  as.data.frame() %>% 
  select(num_beneficiarios_mean, suma_valor, pib)

##      num_beneficiarios_mean  suma_valor          pib
## 0%                   2.0000      3216.0     15515.27
## 25%                346.7083    704285.5     87586.45
## 50%                947.5000   2101369.0    190970.75
## 75%               2439.0417   5734926.8    510100.68
## 100%            451383.8333 853921396.0 763805984.80

¡También podemos revisar los resultados de skimr!

# por grupo:
df %>% 
  group_by(nome_regiao) %>% 
  summarise(
    n = n(),
     expo_pc = mean(exportacoes_pc, na.rm=T)
    ,expo_pc2 = weighted.mean(
      exportacoes_pc, w=populacao, na.rm=T)
    ,popM = sum(populacao, na.rm=T) / 1e6
    ,area = sum(area_geografica, na.rm=T)
  ) 

## # A tibble: 5 x 6
##   nome_regiao      n expo_pc expo_pc2  popM     area
##   <chr>        <int>   <dbl>    <dbl> <dbl>    <dbl>
## 1 Centro-oeste   467   2132.    1601.  16.1 1557857.
## 2 Nordeste      1794    159.     291.  60.6 1630331.
## 3 Norte          450    704.    1214.  16.7 3796198.
## 4 Sudeste       1668    948.    1302.  87.2  919671.
## 5 Sul           1191    669.    1526.  30.5  580495.

# Cuales son las expo. per capita de cada region?

Relaciones entre variables

En general usamos diagramas de dispersión (scatter plots):

ggplot(df, aes(x=suma_perbenef, y=pib_percapita)) +
  geom_point() +
  NULL

ggplot(df, aes(x=log10(suma_perbenef), y=log10(pib_percapita))) +
  geom_point(alpha=0.3, size=0.5) +
  geom_smooth() + # smoother / suavizado
  NULL

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(df, 
       aes(x=log(suma_valor_pc), y=log(pib_percapita), 
           color=semiarido)) +
  geom_point(size=0.5, alpha=0.5) +
  geom_smooth() + # smoother / suavizado
  NULL

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Una relación lineal en un gráfico log-log implica elasticidad constante: la pendiente de un ajuste lineal indica aproximadamente cuánto crece en promedio porcentualmente una variable cuando la otra crece 1%.

# a veces es mejor un plot por grupo:
ggplot(df, aes(x=num_beneficiarios_mean_pc, y=log10(pib_percapita))) +
  geom_point(size=.7, alpha=0.3) +
  geom_smooth(se=F, method="loess", size=.8) +
  facet_wrap(vars(nome_regiao), scales="free_y") +
  theme_light() +
  NULL

## `geom_smooth()` using formula 'y ~ x'

Cuando N es grande, podemos usar alguna de las siguientes opciones:

# partir en cuadrados:
ggplot(df) +
  geom_bin2d(aes(log10(suma_valor_pc), log10(pib_percapita))) +
  scale_fill_viridis_c()

# partir en hexagonos:
# install.packages("hexbin")
ggplot(df) +
  geom_hex(aes(log10(suma_valor_pc), log10(pib_percapita))) +
  scale_fill_viridis_c()

# usar boxplots:
ggplot(df, aes(log10(suma_valor_pc), log10(pib_percapita))) +
  geom_boxplot(aes(group=cut_number(log10(suma_valor_pc), 10)))

  # todos los grupos tienen igual N

# despegar puntitos pegados (con cuidado, p/variables discretas):
# ggplot(df) + 
  # geom_jitter(aes(log10(suma_valor_pc), log10(pib_percapita)))

Correlaciones

Con un correlograma visualizamos la correlación entre todas las variables cuantitativas:

GGally::ggcorr(
  df_num, method=c("pairwise","spearman"),  
  label=T, hjust=1, label_size=2, layout.exp=10, size=3)

Las correlaciones de a pares pueden ser útiles para detectar variables redundantes:

cor_matrix = cor(df_num, method="spearman", use="pairwise")
cor_matrix[upper.tri(cor_matrix, diag=T)] = NA
df_cor = cor_matrix %>% as.table() %>% as.data.frame()
df_cor %>% 
  rename(corr = Freq) %>% 
  filter(!is.na(corr) & Var1 != Var2) %>% 
  arrange(-abs(corr)) %>% 
  head(10) %>% 
  knitr::kable() %>%
  kableExtra::kable_styling()

Var1	Var2	corr
despesa_ciencia_tecnologia_pc	despesa_ciencia_tecnologia	0.9998316
importacoes_pc	importacoes	0.9961378
exportacoes_pc	exportacoes	0.9932097
suma_valor	num_beneficiarios_mean	0.9908110
suma_valor_pc	num_beneficiarios_mean_pc	0.9829813
despesa_pessoal_encargos_sociais	despesa_educacao_cultura	0.9586380
pib	vab_industria	0.9265616
despesa_transporte_pc	despesa_transporte	0.8852189
populacao	despesa_educacao_cultura	0.8561423
populacao	despesa_pessoal_encargos_sociais	0.8551407

Métodos de manejo de NA en cor():

• "everything": usa todas las observaciones (si hay NA en un par de variables, arroja NA)
• "pairwise": usa las observaciones completas de a pares (descarta las observaciones con algún NA de a pares)
• "complete.obs": usa las observaciones completas del data.frame (descarta las observaciones con algún NA considerando todas las filas)

En la clase que viene vamos a ver cómo medir asociaciones entre distintos tipos de variables, y cómo esto nos puede ayudar para guiar el EDA.

Análisis de variables categóricas

Revisamos la distribución de frecuencias de las variables categóricas.

table(df_cat$tipologia)

## 
## IntermediarioAdjacente    IntermediarioRemoto         RuralAdjacente 
##                    686                     60                   3040 
##            RuralRemoto                 Urbano 
##                    323                   1456

table(df_cat$tipologia, df_cat$semiarido)

##                         
##                          FALSE TRUE
##   IntermediarioAdjacente   525  161
##   IntermediarioRemoto       54    6
##   RuralAdjacente          2162  878
##   RuralRemoto              240   83
##   Urbano                  1322  134

# para todas las variables:
tabs = df_cat %>% map(function(x) table(x, useNA="ifany"))
tabs$hierarquia

## x
##   Capital Regional     Centro de Zona       Centro Local Centro Subregional 
##                189                561               4479                164 
##          Metrópole 
##                177

tabs$tipologia

## x
## IntermediarioAdjacente    IntermediarioRemoto         RuralAdjacente 
##                    686                     60                   3040 
##            RuralRemoto                 Urbano                   <NA> 
##                    323                   1456                      5

Podemos hacer tablas más informativas con janitor::tabyl().

df_cat %>% 
  tabyl(tipologia, semiarido) %>% 
  adorn_totals() %>% 
  adorn_percentages("row") %>% 
  adorn_pct_formatting() %>% 
  adorn_ns()

##               tipologia         FALSE         TRUE
##  IntermediarioAdjacente  76.5%  (525) 23.5%  (161)
##     IntermediarioRemoto  90.0%   (54) 10.0%    (6)
##          RuralAdjacente  71.1% (2162) 28.9%  (878)
##             RuralRemoto  74.3%  (240) 25.7%   (83)
##                  Urbano  90.8% (1322)  9.2%  (134)
##                    <NA> 100.0%    (5)  0.0%    (0)
##                   Total  77.3% (4308) 22.7% (1262)

Cortamos el target numérico en bins y hacemos tablas de frecuencias:

tablas_vs_target = function(target, df_cat, variable, g=5) {
  df_cat$target_ = Hmisc::cut2(target, g=g)
  tab = df_cat %>% 
    tabyl(target_, !!sym(variable)) %>% 
    adorn_totals() %>% 
    adorn_percentages("col") %>% 
    adorn_pct_formatting() %>% 
    adorn_ns()
  return(tab)
}

tablas = list()
for (col in names(df_cat)) {
  tablas[[col]] = tablas_vs_target(
    df$pib_percapita, df_cat, col)
}

tablas$hierarquia

##         target_ Capital Regional Centro de Zona  Centro Local
##  [ 4483,  9622)       2.1%   (4)     9.8%  (55)  23.4% (1046)
##  [ 9622, 14545)       5.8%  (11)    19.4% (109)  21.3%  (953)
##  [14545, 22225)      19.0%  (36)    16.4%  (92)  20.0%  (894)
##  [22225, 33591)      25.9%  (49)    27.8% (156)  18.8%  (841)
##  [33591,464883]      47.1%  (89)    26.6% (149)  16.6%  (745)
##           Total     100.0% (189)   100.0% (561) 100.0% (4479)
##  Centro Subregional    Metrópole
##          3.0%   (5)   2.3%   (4)
##         11.6%  (19)  12.4%  (22)
##         24.4%  (40)  29.4%  (52)
##         26.2%  (43)  14.1%  (25)
##         34.8%  (57)  41.8%  (74)
##        100.0% (164) 100.0% (177)

tablas$semiarido

##         target_         FALSE          TRUE
##  [ 4483,  9622)   9.9%  (425)  54.6%  (689)
##  [ 9622, 14545)  16.1%  (695)  33.2%  (419)
##  [14545, 22225)  23.4% (1009)   8.3%  (105)
##  [22225, 33591)  25.2% (1085)   2.3%   (29)
##  [33591,464883]  25.4% (1094)   1.6%   (20)
##           Total 100.0% (4308) 100.0% (1262)

Podemos hacer tablas más lindas con knitr::kable() y kableExtra::kable_styling().

tablas$hierarquia %>% 
  knitr::kable() %>% 
  kableExtra::kable_styling(
    bootstrap_options=c("striped", "hover"), font_size=12
  )

target_	Capital Regional	Centro de Zona	Centro Local	Centro Subregional	Metrópole
[ 4483, 9622)	2.1% (4)	9.8% (55)	23.4% (1046)	3.0% (5)	2.3% (4)
[ 9622, 14545)	5.8% (11)	19.4% (109)	21.3% (953)	11.6% (19)	12.4% (22)
[14545, 22225)	19.0% (36)	16.4% (92)	20.0% (894)	24.4% (40)	29.4% (52)
[22225, 33591)	25.9% (49)	27.8% (156)	18.8% (841)	26.2% (43)	14.1% (25)
[33591,464883]	47.1% (89)	26.6% (149)	16.6% (745)	34.8% (57)	41.8% (74)
Total	100.0% (189)	100.0% (561)	100.0% (4479)	100.0% (164)	100.0% (177)

Y aún más lindas con gt o gtsummary :)

Visualizaciones

En general usamos gráficos de barras:

ggplot(df) +
  geom_bar(aes(x=nome_regiao)) +
  NULL

Vamos a mejorarlo un poco:

# ordenando por N
ggplot(df) +
  geom_bar(aes(x=fct_infreq(nome_regiao))) +
  NULL

Agregando labels:

# La version mas general: DF agrupado + geom_col()
# (en lugar de geom_bar())
df_tmp = df %>%
  group_by(nome_regiao) %>% 
  summarise(
    n = n()
  ) %>% 
  mutate(
    porcentaje = n / sum(n) * 100,
    nome_regiao = fct_reorder(nome_regiao, n)
    )
ggplot(df_tmp) +
  geom_col(aes(x=nome_regiao, y=porcentaje), alpha=0.5) +
  geom_text(aes(
    x=nome_regiao, y=porcentaje, label=glue("{n}\n({round(porcentaje,2)}%)")), 
            color="black", size=3) +
  theme_light() +
  labs(x="Region", y="N") +
  NULL

Nota: a veces con una buena tabla alcanza!

Podemos analizar la relación con alguna variable numérica usando alguna medida resumen como la media:

ggplot(df) +
  geom_bar(aes(x=nome_regiao, y=pib_percapita), 
           fun="mean", stat="summary", alpha=.5) + 
  theme_minimal() + 
  NULL

Los puntos suelen usarse para representar valores donde no es necesario incluir el 0.

ggplot(df, aes(x=nome_regiao, y=pib_percapita)) +
  geom_point(fun="mean", stat="summary") +
  # stat_summary(fun = "mean", geom = "point") +
  geom_line(aes(group=1), fun="mean", stat="summary") +
  theme_minimal() + 
  NULL

Ordenando los niveles según los valores de otra variable:

# resumir los datos y luego usar geom_col():
df_tmp = df %>% 
  group_by(nome_uf) %>% 
  summarise_at(vars(populacao), sum)
ggplot(df_tmp) +
  geom_col(
    aes(x=populacao, y=fct_reorder(as.factor(nome_uf), populacao)), 
    alpha=0.5) +
  labs(y=NULL) +
  scale_x_continuous(labels=scales::comma) +
  theme_light() +
  NULL

# # La otra manera:
# ggplot(df) +
#   geom_bar(
#     aes(x=populacao, y=fct_reorder(as.factor(nome_uf), populacao, .fun=sum)), 
#     fun="sum", stat="summary", alpha=0.5) +
#   labs(y=NULL) +
#   scale_x_continuous(labels = scales::comma) +
#   theme_classic() +
#   NULL

También podemos hacer un boxplot por grupo, como vimos antes.

Para pensar

Vamos a suponer momentáneamente que nuestro target es binario: PIB bajo (0) o PIB alto (1). ¿Qué tablas y visualizaciones nos pueden servir como EDA en este escenario?

df = df %>% 
  mutate(target_pib = ifelse(pib_percapita > mean(pib_percapita), 1, 0)) 

# TODO

Conclusiones

• El EDA suele ser un proceso iterativo (explorar → transformar → explorar) que fácilmente se puede volver caótico. Por eso es importante:
  • Tener un plan más o menos claro de los pasos que vamos a hacer, definiendo a qué información indispensable queremos llegar
  • Anotar durante el proceso los nuevos análisis e hipótesis que se desprenden
  • Documentar todo en la notebook/rmarkdown/script
  • Seleccionar los cuadros y gráficos relevantes para comunicar los resultados, y contar lo importante del proceso
• Algunos aspectos importantes a tener en cuenta antes/durante el EDA:
  • Conocer la fuente de los datos y posibilidad de actualizarlos/replicarlos en el futuro
  • Saber la definición de las variables
  • Revisar registros duplicados / información innecesaria
  • Ser conscientes de qué observaciones eliminamos/conservamos en cada paso
  • En general es aconsejable primero “limpiar los datos” (filtrar observaciones/variables, analizar/imputar valores faltantes, consolidar tablas, etc.) antes de crear nuevos features
  • Conocer la distribución de las variables categóricas/numéricas con estadísticos, tablas y visualizaciones, usando la escala apropiada
  • Conocer cómo se relacionan las variables entre sí, usando métricas y visualizaciones apropiadas
  • Poner especial atención en la distribución del target y su relación con el resto de las variables
  • ¿Lo que observamos es razonable? ¿Hay cosas inesperadas? ¿Hay variables redundantes?
• El objetivo del EDA es entender los datos. Por lo tanto, los pasos de un EDA no necesariamente vayan a ser los mismos que apliquemos cuando entrenemos un modelo supervisado.

Referencias

• R for Data Science
• R Markdown: The Definitive Guide
• RMarkdown for Scientists
• Python Data Science Handbook
• Hands-on Exploratory Data Analysis With Python
• Librerías para Python:
  • ydata-profiling
  • sweetviz

Bonus tracks

Con plotly::ggplotly() podemos hacer interactivo cualquier gráfico de ggplot2.

g = ggplot(
  df
  ,aes(x=log(pib_percapita+1), y=log(exportacoes+1)
       ,color=semiarido, label=municipio, label2=nome_regiao)) +
  geom_point(alpha=0.8, size=0.8) +
  theme_minimal() + 
  NULL
plotly::ggplotly(g)

Con GGally::ggpairs() podemos analizar relaciones entre variables de múltiples tipos.

df_reduced = df %>%
  select(
    nome_regiao, suma_perbenef, densidad,
    despesa_ciencia_tecnologia_pc, pib_percapita, semiarido
  ) %>% 
  mutate_if(is.numeric, function(x) log10(x+1))
GGally::ggpairs(
  df_reduced, 
  columns = c("suma_perbenef", "densidad",
            "despesa_ciencia_tecnologia_pc", "pib_percapita"),
  mapping = aes(color=semiarido),
  lower = list(continuous = GGally::wrap("points", alpha=0.3, size=0.1))
)

¿Cómo generamos todos los gráficos bivariados con respecto a la variable target?

# generamos log de todas las vars
df_num_expanded = df_num %>%
  mutate_all(list("log" = function(x) ifelse(x==0, 0, sign(x)*log10(abs(x))) ))

# Medimos corr en logs vs el target (equivale a spearman!)
corr_mat = cor(
  df_num_expanded %>% select(ends_with("_log")) %>% select(-pib_percapita_log),
  df_num_expanded %>% select(pib_percapita_log),
  method="pearson", use="pairwise")
df_corr = corr_mat %>% as.data.frame()
names(df_corr) = "correlation"
df_corr = df_corr %>% 
  arrange(-abs(correlation)) %>% 
  rownames_to_column("variable")

# plot function
plt_corr = function(df, variable) {
  p = ggplot(df, aes(x=!!sym(variable), y=pib_percapita_log)) +
    geom_jitter(size=.3, alpha=.2) +
    geom_smooth(se=T, method="loess", size=.8, formula=y~x) +
    labs(title=variable) +
    theme_minimal() +
    NULL
  return(p)
}

# lista de plots
variables = df_corr$variable
list_plots = map(variables, function(x) plt_corr(df_num_expanded, x))

# save plots
ggsave("output/bivariado_pib.pdf",width=6, height=6,
       gridExtra::marrangeGrob(grobs=list_plots, nrow=2, ncol=1))

Mapas

Descargamos los shapefiles de IBGE.

Shapefile es un formato de archivo con información geográfica. Indican el sistemas de coordenadas de referencia y proyección cartográfica usados. Contiene puntos (posiciones), líneas (recorridos) y polígonos (superficies) de alguna parte del mundo.

Los podemos leer con sf::read_sf() y hacer un join con nuestros datos.

sf_mun = sf::read_sf(
  "data/raw/ibge/br_municipios_20200807/BR_Municipios_2019.shp")
sf_uf = sf::read_sf(
  "data/raw/ibge/br_unidades_da_federacao/BR_UF_2019.shp")

names(sf_mun)

## [1] "CD_MUN"   "NM_MUN"   "SIGLA_UF" "AREA_KM2" "geometry"

dim(sf_mun)

## [1] 5572    5

sf_mun = sf_mun %>%
  left_join(df, by=c("CD_MUN"="codigo_municipio"))

# municipios faltantes en nuestros datos
sf_mun %>% 
  filter(is.na(municipio)) %>% 
  select(CD_MUN, NM_MUN, SIGLA_UF, AREA_KM2)

## Simple feature collection with 2 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -53.53194 ymin: -33.61715 xmax: -50.53778 ymax: -30.00031
## Geodetic CRS:  SIRGAS 2000
## # A tibble: 2 x 5
##   CD_MUN  NM_MUN          SIGLA_UF AREA_KM2                             geometry
##   <chr>   <chr>           <chr>       <dbl>                   <MULTIPOLYGON [°]>
## 1 4300001 Lagoa Mirim     RS          2872. (((-53.52638 -33.61715, -53.5284 -3~
## 2 4300002 Lagoa dos Patos RS         10197. (((-52.09772 -32.16174, -52.09888 -~

Creamos mapas coropléticos con geom_sf(). En estos mapas las áreas se sombrean de distintos colores según alguna característica.

sf_mun_tmp = sf_mun %>% filter(nome_uf == "Sergipe")
sf_uf_tmp = sf_uf %>% filter(NM_UF == "Sergipe")

ggplot() +
  geom_sf(data=sf_mun_tmp, aes(fill=pib_percapita), color="black", size=0.5) +
  geom_sf(data=sf_uf_tmp, fill=NA, color="black", size=0.8) +
  scale_fill_viridis_c(trans="log10", na.value="gray60") +
  # coord_sf(xlim=c(-75,-55), ylim=c(-10, 3)) +
  labs(fill="PIB per capita") +
  NULL

g = ggplot() +
  geom_sf(data=sf_mun, aes(fill=pib_percapita), color=NA) +
  geom_sf(data=sf_uf, fill=NA, color="black", size=0.5) +
  scale_fill_viridis_c(trans="log10", na.value="gray60") +
  theme_minimal() +
  guides(scale ="none") +
  theme(legend.position=c(0.15,0.2)) +
  labs(fill="PIB per capita") +
  NULL

ggsave("output/mapa_pibpercapita.png", g, width=10, height=8)

Si tuviéramos lat y long en nuestros datos podríamos inferir el mapa de esta manera:

df_geo = df %>% 
  st_as_sf(coords=c("longitud", "latitud"), crs=4326)

ggplot() + geom_sf(df_geo)